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Abstract 

A recent proposal^ for practical calculation of vibrational mode lifetimes is tested on simple, 
low-dimensional anharmonic models. The proposed scheme approximates the mode lifetime in 
terms of ensemble averages of specific functions in phase-space; various levels of approximation 
correspond to ensemble moments of the Liouvillian. It is shown that, for systems where the 

vibrational density of states is well- approximated by a single broadened peak, the fourth-moment 

£ 
i approximation works well over the full range of temperature. 
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I. INTRODUCTION 

Dickel & Daw^ recently proposed an efficient, approximate means of calculating vi- 
brational mode lifetimes in solids. The method involves ensemble averages of appropriate 
functions in phase space that can be carried out by conventional Monte Carlo in combi- 
nation with a means of calculating forces, such as interatomic potentials or first-principles 
electronic structure codes. The approach was illustrated on a lattice model of non-linear 
interactions, where the dependence of the mode lifetimes on cell size and temperature was 
investigated numerically. 

While the aim of the original work was to further calculations of vibrational mode lifetimes 
in solids, the purpose of the present work is to examine in more detail the approximations 
involved in the method. To this end we take up the same method as applied to very simple 
systems of just one or two degrees of freedom. In considering systems of such simplicity, we 
analyzed some aspects of the problem analytically as well as numerically, and the insights 
obtained are reported here. These insights are expected to prove fruitful in the application 
of this method to the original target (vibrational lifetimes in solids). 

This paper is organized as follows. First, we recap briefly the proposal of Dickel & Daw 
(DD). Then we apply the proposed method to the simple dynamical models considered here. 
Our analysis of the results focuses on the density of states, by which we can understand when 
and why the approximations work as they do. Finally, we draw our conclusions. 

II. BACKGROUND & SCOPE OF THE PRESENT WORK 

We summarize here the proposed method of DD, who began by examining the momentum 
Auto-Correlation Function (MACF) 

where the angular brackets indicate phase-space averages over the canonical ensemble at 
temperature T (p — exp (-H/T)). 

The auto-correlation can be studied in terms of the LiouvilliarP^, which governs the time 
evolution of functions f(x,p,t) in phase space according to 

df . t . 



where the (hermitian) Liouvillian operator is 

,dH d dH d 



L = t{H,} = tJ2( 



dxi dpi dpi dxi 
The equation of motion can be integrated formally, so that 

f(x,p,t) = e- ul f(x, P ,0) 

and we can express the auto-correlation explicitly in terms of L: 

_ (pe~ itL p) 
X[) ~ (P 2 ) 

The Taylor Series of x(t) 

t 2 t 4 t 6 

X{t) = l ~^2\ +/i4 4! _/i6 6! + '" 

relates the derivatives of x{t) a t t = to the moments of the Liouvillian acting on the 

momentum: 

(pL n p) 

These moments are also the moments of the density of states (DOS) derived from x(t). That 
is, taking the Fourier transform of x{t) to get n(u), the moments are also 

/+oo 
duj u m n(u) 
-oo 

Auto-correlation functions considered in this work typically have strong oscillations damp- 
ened by some sort of decaying envelope (for examples, see Figs, lip])- We propose here to 
use the area under the square of the MACF as a measure of the lifetime 

"+0O 



dt X (t) 2 (2) 

This is not intended to correspond to a particular physical measurement that might be 
performed, but rather is suggested as a simple generic measure of the rate of the decay of 
the correlation. Such a measure also lends itself easily to analysis, as we shall see. Using 
Parseval's Theorem, the lifetime is also given as the area under the n{oj) 2 curve: 



r = / dt x(t) 2 = dco n(u) 2 (3) 
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DD observed that the lifetime r can be expressed as a function of the moments 

T = F(/i 2 ,/i 4 ,/l 6 ,...) 

which can be re-expressed (using dimensional analysis) as 

t/t 2 = G(7 4 , 7 6 , . . .) 

where r 2 = //J and the 7's are dimensionless parameters 



(/i 2 )"/ 2 

that characterize the shape of the DOS. While it is not generally possible to know all of 
the moments, DD proposed that in certain circumstances the lifetime might be practically 
approximated from a knowledge of only the lowest moments. This suggests a series of 
approximations, starting with only the second moment 

r = ct 2 (4) 

and including successively higher moments. The fourth moment approximation would then 
be 

r = r 2 F( 74 ) (5) 

where F is some function yet to be determined. The higher moments correspond to ensemble 
averages of higher powers of the Liouvillian, and so each higher moment involves higher time 
derivatives of the dynamical variables. 

DD then went on (in part 2) to test the lowest approximation on a simple model of 
non-linear lattice vibrations as a function of cell size and temperature. First, much as 
done by Ladd, et a/.pDD calculated from ordinary molecular dynamics the auto-correlation 
function for each normal mode a periodic cell of a given size (appropriately sampled at the 
specified temperature) and from there the lifetime. Second, they calculated using standard 
Monte Carlo the second moment /z 2 (hence r 2 ) for each mode. (This second part of the 
demonstration is, of course, requires much less computational time than the first.) They 
then plotted t/t 2 vs. temperature for all modes, and found that at high temperatures 
the lifetime was simply proportional to r 2 . Furthermore, at high temperature, the auto- 
correlation functions scaled in a simple way. That is, plotting all of the calculated x vs - 
t/r% exhibited a data collapse, revealing that indeed the high-temperature dynamics of the 



mode decay could simply be described by a single parameter. Thus, the high temperature 
behavior was well approximated at the lowest level (second moment). 

DD ended by speculating that the behavior over the full range of temperature might be 
accounted for by including fourth moment, but that was not tested. Also, that paper did 
not offer much insight as to why the second moment approximation should work well at high 
temperature but be insufficient at low temperatures. 

The present study uses several simple dynamical models as the basis for testing the fourth 
moment approximation and also in using the density of states to provide an analysis of why 
the approximation might work and when it would be expected to fail. 

III. MODELS CONSIDERED 

We consider three simple model hamiltonians in one (x) and two (x and y) dimensions. 
These models are chosen because they are simple, non-linear, and the ensemble averages 
can be obtained analytically. The momentum conjugate to x is p\ that to y is q. 

x 4 model: 

H(p, x) = p 2 + x 2 + x 4 (6) 

The auto-correlation in the x 4 model has been studied extensively before™. In that work, 
an analytic approximation to x{t) was obtained at low temperature: 

_ cos(t) -3Ttsin(t) 

showing an oscillatory behavior with an algebraically decaying envelope. Our calculated 
auto-correlation conforms well to this analytical form at low temperatures. 

x 2 y 2 model: 

H(p, x, q, y) =p 2 + q 2 /M + x 2 + y 2 + x 2 y 2 (8) 

The x 2 y 2 model is a simple extension to two modes coupled nonlinearly In this model, 
we investigate various values of the ratio (M) of the masses between the two modes, which 
controls the degree of resonance. 



"cubic" model: 



H(p, x, q, y)=p 2 + q 2 + x 2 + y 2 + -{x 2 + y 2 ) 2 + -(x 3 



3xy 2 



(9) 



The "cubic" model for certain parameters has multiple minima in the x — y plane and 
exhibits a "structural" transformation (from multiple attractors to a single attractor) with 
temperature, which makes it interesting to include in the present study. To explore the 
effects produced by this transition, we tried various values of the strength (A) of the cubic 
term. For A < 2/9, there are 3 off-center global minima with one local minimum on-center. 
For A > 1/4 only there is only 1 global minimum on-center. 



Some examples of a calculated MACF are shown in Figs. [TJJ2] For the x 4 model, the 
function exhibits a simple oscillation and decay. In the "cubic model" , the function displays 
less regularity because of the less symmetrical potential. 
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FIG. 1. The MACF at three temperatures for the x 4 model. 



o 




FIG. 2. The MACF of the x-mode at A = 0.5 and T = 0.2 for the "cubic" model. 



IV. TESTING THE FOURTH MOMENT APPROXIMATION 



We want to determine if the form in Eq. [5] is robust enough to approximate the lifetimes 
in the various simple models we have chosen. In the x A model, for example, we can perform 
ensemble dynamics at various temperatures and extract the lifetime by Eq. [5j The lifetime 
vs. temperature for the a; 4 model is then shown in Fig. |3J 




FIG. 3. Lifetime (Eq. 0) vs. temperature in the x 4 model. 



In view of Eq. M we represent these results as a scatterplot of rjr^ vs. 74, where the 
temperature-dependent r 2 and 74 are calculated analytically. Noting that 74 > 1, and the 
power-law behavior of r and the moments with T, we will plot Xogir/T^) vs. log(74 — 1). 
This is shown in Fig. [4] 
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FIG. 4. Scatterplot of t/t2 vs 74 for the x model. The straight line is a fit using Eq. 10 



Similar results can be seen for the x y model (see Fig. 
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FIG. 5. T/V2 vs 74 of the x-mode for the x y model. The straight line is a fit using Eq. 10 



The relations in Figs. |4]j5] are fit well by 

r/r 2 ex (74 - ir 1/2 (10) 



From the analytic approximation to MACF for the x 4 model at low temperature (Eq. u\ 
we can also calculate r, /x 2 , and /x 4 , and we find 



/r 2 = tt( 74 - l)- 1 / 2 (11) 



T 



The analytical form for x(t) was derived by Sen, et al. only for low temperatures, by noting 
the dependence of the oscillator frequency on energy and the contributions of different 
energies in the canonical ensemble. However, here we find the relationship between r, r 2 , 
and 74 extends over a large range of temperature. The reason for this extended range will 
be understood better below. 

At low T, the low moments for both x 4 and x 2 y 2 models behave similarly, in that yU 2 ~ 
1 + aT and /z 4 ~ 1 + 1aT so that 7 4 approaches 1 as T 2 . Thus r 2 is approaches a constant 
while 74 — 1 goes to zero, and the lifetime diverges like r « T~ 4 at low temperature. The 
temperature dependence at low T is dominated by the approach of 74 to 1. 

At high temperature, the moments for the x 4 model go as /i 2 ~ aT 1 / 2 + b and /z 4 ~ 
cT — dT 1 ! 2 . So 7 4 saturates to a constant as T _1//2 , leaving only the variation in r 2 to 
account for the change in lifetime. Thus the lifetime at high T, is governed by the behavior 
of r 2 and r « T~ 1 / 4 . 

This accounts well for the two power-law regimes visible in Fig. [3j 

For the x 2 y 2 model, by contrast, at high temperature, the moments go as /i 2 ~ 2T / logT+ 
1/2 and /i 4 ~ AT 2 / log T+4T/M which makes 74 go as log T. This cancels a log T dependence 
in r 2 leaving r rj T -1 / 2 . 

The x 4 and x 2 y 2 models seem to be well-described by the simple combination of the first 
two moments. However, by contrast, the corresponding scatterplot for the "cubic model" 
deviates significantly (Fig. |6|), so that there is no simple relationship between r and the 
first two moments. Evidently, higher moments will be required to capture the dynamical 
behavior of the cubic model over a wide range of temperatures and parameters. 
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FIG. 6. Scatterplot of T/V2 vs 74 of the x-mode for the "cubic" model, showing irregular behavior 
as compared to the other models (Figs. (4l5|. 



V. ANALYSIS 



From the previous results, we see that the behavior of the lifetime for the x 4 and x 2 y 2 
models over a wide range of parameters and temperature is captured in the behavior of the 
two lowest moments (^2 and /X4) which can be calculated analytically. However, for the cubic 
model, the behavior is more complex, requiring at least higher moments in the description. 
We investigate here the reasons for success in one case and not in the other. 

Fig. u\ shows the insight gained from checking for a data-collapse for x(^)> by scaling the 
time t by the lifetime r (Fig. pi) for the x 4 model. The results illustrate that while the 
oscillations of auto-correlation functions vary with temperature, they are contained by one 
decaying envelope, which is what we are trying to capture. 
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FIG. 7. Data collapse of MACF for the x 4 model (as explained in the text.) 



As one might expect from the data-collapse, the DOS for the x 4 model is also simple, as 
shown in Fig [8] for various temperatures. 
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FIG. 8. Density of states at various temperatures for the x 4 model. 



The DOS of this model is characterized by a single dominant peak that shifts and broadens 
with temperature, as one would typically expect of a vibrational mode in an anharmonic 
solid. In such a case, the lifetime depends mostly on the shape of the DOS around the peak, 
and two parameters (peak value of the DOS and the width) are sufficient to describe it. At 
low temperatures, 74 — > 1, while at high temperatures 74 — > 2.2 (for this model). Recalling 
74 as the (dimensionless) ratio 114/ 'ji 2 ,, it is aptly designated as a "shape parameter" of the 
DOS. 

The DOS of the x 2 y 2 model (Fig. |M) is only somewhat more complex than that of the x 4 
model. 
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FIG. 9. Density of states of the y-mode at M = 1 for the x y model 

The simple evolution of the DOS with the temperature and other parameters for these 
models explains why a simple, generic relationship can exist between r and the first two 
moments of the DOS. To explore this point further, let us consider a generic, single-mode 
DOS that is peaked at an oscillator frequency uq and broadened to a width Q. Both the 
oscillator frequency and width will depend on temperature. At low temperatures, Q « ojq, 
and from Eq. [3] we have 

The leading behavior of the lowest two moments is 

/i 2 « wjftl + an 2 J 'ul) 

/i 4 w w£(l + btt 2 /ul) 
where a and b depend on the details of the DOS. Then 



and 



T/V 2 PS Uq/VL 



7 4 - 1 « n 2 /ul 
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Eliminating Q and loq among the two relations gives 

~/r 2 » ( 74 - I)- 1 ' 2 



T I 



just as we found in Eq. 10 So long as the DOS has this simple, generic behavior, the same 
relationship obtained here should hold. 

At high temperatures, if the DOS can be assumed to be a mostly featureless and broad 



distribution with width Q and height f2 _1 , then r 
the shape parameter saturates at some value (74 f 
tracked by that of r 2 , so that 

T ^ T 2 

which is the behavior reported by DD. 



: Q x and /i 2 ~ ^ 2 so r 2 « Q 1 . While 
c), in which case the variation in r is 



The DOS of cubic model (Fig. 10) is much more complicated than that of x 4 and x 2 y 2 
model, which explains why the simple 2-parameter scatterplot (Fig. [6]) does not capture the 
behavior. 
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FIG. 10. Density of states of the x-mode at A = 0.2 and various temperatures for the cubic model. 

Finally we note that 74, in addition to being a simple measure of the shape of the DOS, 
is also a direct measure of the degree of anharmonicity of the mode as averaged over the 
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ensemble. Specifically, 74 for a given mode can be re-written as 

71 = -w (12) 

where / is the force associated with a displacement x. A harmonic system is, of course, 
defined where the force obeys / + kx — 0. In the anharmonic ensemble, we could define an 
effective k by that which minimizes the deviation from linear. That is, define the effective 
k by minimizing a = ((/ + kx) 2 ). The minimum value of a then measures the degree of 
anharmonicity of the system as effective for the ensemble. For a harmonic system, a m ; n = 0. 
In general, k e & = —(xf)/(x 2 ) and 

«min = -r^rili ~ 1) ( 13 ) 

[X ) 

showing how the deviation 74 — 1 is directly related to the effective anharmonicity of the 
ensemble. 



VI. CONCLUSIONS 

We have investigated using low- dimensional models the proposal that the mode lifetime 
in equilibrium might be approximated from the two lowest moments of the Liouvillian. For 
the generic case of a DOS dominated by a single peak broadened and shifted, as is the case 
here for the x A and x 2 y 2 models, we see that the fourth moment approximation works well. 
In the case of the cubic model, the fourth moment approximation is insufficient, which can 
be understood in terms of the more complex structure of the DOS. The multiple minima of 
the cubic model creates a more complex dynamics that cannot be captured with just two 
parameters. 
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